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The triumph of effective mass theory in describing the energy spectrum of dopants does not guar¬ 
antee that the model wavefunctions will withstand an experimental test. Such wavefunctions have 
recently been probed by scanning tunneling spectroscopy, revealing localized patterns of resonantly 
enhanced tunneling currents. We show that the shape of the conducting splotches resemble a cut 
through Kohn-Luttinger (KL) hydrogenic envelopes, which modulate the interfering Bloch states of 
conduction electrons. All the non-monotonic features of the current profile are consistent with the 
charge density fluctuations observed between successive {001} atomic planes, including a counter¬ 
intuitive reduction of the symmetry - a heritage of the lowered point group symmetry at these 
planes. A model-independent analysis of the diffraction figure constrains the value of the electron 
wavevector to ko = (0.82 ± 0.03)(27r/asi). Unlike prior measurements, averaged over a sizeable den¬ 
sity of electrons, this estimate is obtained directly from isolated electrons. We further investigate 
the model-specific anisotropy of the wave function envelope, related to the effective mass anisotropy. 

This anisotropy appears in the KL variational wave function envelope as the ratio between Bohr 
radii h/a. We demonstrate that the central cell corrected estimates for this ratio are encouragingly 
accurate, leading to the conclusion that the KL theory is a valid model not only for energies but for 
wavefunctions as well. 


I. INTRODUCTION 

Modern applications of electronic quantum control at 
the single dopant level underline significant scientific 
challenges to the theory of doped semiconductors. Many 
challenges revolve around the determination of the donor 
electron wavefunction. In silicon, the rich conduction 
band structure renders the problem of hydrogenic impu¬ 
rities unsolvable even within the simple effective mass 
approximation. The pioneering work of Kohn and Lut- 
tinger (KL) provides a tentative answer within a vari¬ 
ational framework [1]. It is unclear, though, if the sim¬ 
plicity of effective mass approximation can withstand the 
comparison with the directly probed charge distribution 
of a donor. 

Furthermore, the variational principle guarantees that 
the minimal energy mean value is closest to the real 
ground state energy, but ascertains nothing about the 
optimal variational wavefunction. Variations with re¬ 
spect to the true ground state wavefunction only lead 
to second-order shifts in the ground state energy, mean¬ 
ing that even sizeable deviations from the ground state 
wavefunction can still lead to reasonably accurate en¬ 
ergies. For example, energies estimated from the KL 
trial function with two parameters are only less than 1% 
higher than a model with four parameters by Kittel and 
Mitchell [2]. 

This indeterminacy represents a threat to the set of 
ideas that have led the pursuit of quantum technolo¬ 
gies. To tame the electron is to manipulate its wave- 
function. Entanglement between electrons is thought to 
be achievable by simply overlapping their wavefunctions. 


as to explore the Pauli exclusion principle and the spin- 
spin effective interaction that comes with it. Knowledge 
largely based on the effective mass KL model for dopant 
wavefunctions is the common root supporting all designs 
of donor-based quantum devices in Si [3]. Should this 
wavefunction be sizeably wrong, the feasibility of most 
quantum devices would need to be revisited. 

Recent scanning tunneling microscopy (STM) images 
of single donors near a silicon (001) surface revealed the 
theoretically predicted intricate interference patterns of 
such electronic states. [4] At first sight, these images do 
not communicate the simplicity of the KL wavefunction 
- two kinds of non-trivial images are revealed, which we 
call here butterfly (B) or caterpillar (C). In Ref. [4], these 
are referred to as types A and B, respectively. The B- 
type has a nodal line akin to a p orbital symmetry, while 
the C-type has a symmetry closer to s. Both patterns are 
diagonally aligned and present mirror symmetry with re¬ 
spect to either [110] or [110] surface directions. The sym¬ 
metry lines are consistently orthogonal to the dimer rows 
direction on the surface (which run along either direction 
in different terraces). 

Here we reconcile the intrincate current profiles of the 
STM image with the simple KL theory by carefully con¬ 
sidering the crystalline and electronic structures of bulk 
silicon. Perturbations by the STM tip, surface relax¬ 
ation, reconstruction or passivation are not considered. 
The model is briefly reviewed in Sec.|n| and the notation 
for the atomic planes is discussed. Section [Hi] juxtaposes 
theoretical and STM images for donors buried at differ¬ 
ent depths. In Sec. El we theoretically dissect the roles 
of the anisotropic mass, Bloch functions, and periodic¬ 
ity of the lattice, and thus pinpoint the main ingredients 
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driving the anisotropy of the observed images. The ori¬ 
gin of the low-frequency interference patterns is detected 
in Sec. |V| analysing the charge distribution in Fourier 
space. This /c-space analysis is further explored in order 
to estimate and h/a from the donor bound states in 
Sec. ED Section Ivnl is devoted to our conclusions and 
final remarks. 


II. MODEL 

Our model is formulated in the context of well estab¬ 
lished bulk Si structural and electronic properties. Sili¬ 
con crystallizes in the diamond structure which consists 
of two interpenetrating face centered cubic (FCC) lat¬ 
tices shifted by 1/4 of the cubic conventional cell body 
diagonal. All lengths referring to atomic positions and 
distances are given here in units of the cubic cell lattice 
parameter asi = 0.543 nm. 

Fig. 0 a) shows a unitary FCC cube where atomic po¬ 
sitions of the diamond structure are indicated. We iden¬ 
tify the stacking of 4 inequivalent (001) atomic planes at 
heights 2 ; = 0,1/4,1/2, 3/4 which we call Ao,Ai/ 4 ,Ai /2 
and A3/4 or, in general Aj/4 with j = 0,1,2,3. In the 
perfect crystal, the next plane above, at z = 1 (also in 
the figure) is equivalent to the Aq plane. It is convenient 
to label planes intercalated midway between consecutive 
A—planes. We call Ij /4 the intercalated plane a distance 
1/8 above Aj/ 4 . 

In the event that a substitutional donor (which defines 
the coordinates origin) is located a distance d from a 
( 001 ) surface at 2 ; = d, translational symmetry is lost. 
We then refer to a general atomic plane at 2 : = n + jf/4 
as with the donor at the Aq^^ plane. Consistently, 
intercalated planes at 2 ; = n + jf/4 + l /8 are called . 

Nearest neighbor (consecutive) atomic planes belong 
to different FCC sublattices. The tetrahedral bonds for 
atom pairs (one in each plane) define zig-zag paths cross¬ 
ing a single /—plane, as illustrated in Fig. Bb). Note that 
bonds across /j^] for j =0 or 2 are overall aligned with the 
[110] diagonal while for j = 1 or 3 the paths follow the 
[110] diagonal. This applies to all j and (n) planes, in¬ 
cluding the surface at 2 ; = d, thus giving rise 

to two possible directions for the atomic reconstruction 
(dimerisation) on a Si (001) surface. In our theoreti¬ 
cal analysis, we infer the surface dimer rows orientation, 
which is the same as the zig-zag overall direction, ac¬ 
cording to the value of j = 4(d — n) with n integer and 
j = 0, C 2, 3. 

Fig. Ifb) also illustrates that the symmetry of the A^y^ 
plane depends on the sublattice to which this plane be¬ 
longs. For j = 0, 2 - indices related to the same sublat¬ 
tice as the donor - the planes show fourfold symmetry 
while those at the other sublattice j = 1,3 have a lower, 
twofold symmetry. Note that the z-axis crosses an atomic 
site only if it belongs to a j = 0 atomic plane. 


Electronic properties of Si are characterized by a band 
structure with sixfold degenerate conduction band min¬ 
ima (valleys), located at along the (100) directions, 
/i = ± 2 :. The valleys are anisotropic, z.e., 

they present different longitudinal and transverse effec¬ 
tive masses (my = 0.9163me and m± = 0.1905me re¬ 
spectively 0 )- Early measurements [6] estimated the 
wavevector of the minima as |k^| = ko = (0.85 ± 




— ). Our measurements allow a less model- 


0 

dependent, more accurate estimate of /cq, see Sec. [Vr 


The presence of a substitutional donor breaks the 
translational symmetry. In terms of electronic structure, 
a simple and successful description of shallow donors in 
Si was presented by Kohn and Luttinger within effec¬ 
tive mass theory (EMT) [T]. The singular donor poten¬ 
tial couples different valleys, leading to a non-degenerate 
ground state which involves a symmetric combination of 
the six valleys. The donor ground state variational wave- 
function proposed by KL has the correct Ai-symmetry 
and is written in terms of envelopes and Bloch functions 
for each conduction band minimum 


= ( 1 ) 

where F^{t) are envelope functions and u^{t) are the 
periodic parts of the Bloch functions, given explicitly in 
Ref. [7] as obtained within first principles density func¬ 
tional theory. As a consequence of the mass anisotropy, 
the envelope functions have the shape of a deformed Is 
orbital 
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and similarly for the x and y valleys. 

The effective Bohr radii a and b were calculated vari- 
ationally by KL for a Coulomb donor potential. Ad¬ 
ditionally we take into account a central cell correction 
potential which is donor species dependent and chosen to 
reproduce each experimental ground state energy. [8] No¬ 
tice that the phenomenological central cell has a different 
radius Vcc for the anisotropic model of the mass compared 
to the spherical mass model adopted in Ref. [8]. Eor the 
anisotropic model, the radii are rcc(P) = 120 pm and 
Tec (As) = 128 pm for P donors and As donors, respec¬ 
tively. 

Eor P donors we get variational radii a = 1.13 nm and 
b = 0.60 nm, while for As donors we get a = 0.86 nm 
and 6 = 0.46 nm. Consequently, the envelope anisotropy 
is b/a = 0.53 irrespective of the donor species. An 
isotropic approximation for the envelope exp(—r/a) leads 
to a single average Bohr radius a = 1.106 nm for P and 
a = 0.815 nm for As. [8] The exact value of a and b have 
an impact only on the size of the image, not on the de¬ 
tails of its symmetry and oscillations. We adopt from 
now on a = 0.90 nm and b = 0.52 nm for the theoreti¬ 
cal calculations (which lead to the same b/a ratio as the 
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FIG. 1. (color online) Bulk crystal structure and notation adopted for atomic and interstitial (001) planes, (a) Diamond 
structure sites represented by dots; an FCC cubic unit cell is also shown. The four inequivalent atomic planes ^o, ^i/ 4 , ^ 1 / 2 , 
and As/ 4 , are indicated. Each interstitial plane is labelled as the atomic plane immediately below it, that is, /o, /1/4, /1/2? and 
G/ 4 - (b) Expanded view of the plane showing nearest-neighbors bonds and interstitial planes across them. The red 

bonds (M-shaped) and the blue ones (W-shaped) form zig-zag paths. The donor site and the symmetry axis z (see text) are 
indicated, (c) Overall view of the crystal structure with a substitutional donor impurity located at the Aq^^ plane. Each cube 
here is a replica of the one in frame (a). (d) Details of the charge distributions of a combination of Bloch functions pinned at 
the donor site. Cuts are shown at the interstitial planes and One may identify shapes anticipating the B and C forms 
observed experimentally. 


original KL theory but incorporates the radii reduced by 
the central cell). 

Our simulated STM images are generated using the fol¬ 
lowing conditions: (i) Tersoff-Hamann approximation, [9] 
in which the tunneling current is proportional to the lo¬ 
cal density of states integrated over the bias energy win¬ 
dow; (ii) constant-height mode; and (iii) the assumption 
that the electronic density associated to the defect level 
is not substantially modified by the presence of the sur¬ 
face, meaning also that the effects of surface reconstruc¬ 
tion and hydrogen saturation are disregarded. Assump¬ 
tion (iii) seems rather drastic, and its validity can only 
be assessed a posteriori, by direct comparison to the ex¬ 
perimental STM images. A justification to attempt it 
comes from the agreement between tight-binding results 
and experiments in Ref. [4], where the former indicates 


that valley populations for donors > 2.5 nm from the 
surface differ from the bulk values by less than 5%. We 
choose the STM tip height to correspond to a distance 6 
above the surface atomic layer 

Pd{x,y) = \'^{x,y,d + 5)\‘^, (3) 

with S = 1/8, i.e. halfway between atomic planes and 
coincident with the intercalated planes Physically, 

this is justified by the fact that, in our bulk-truncated 
model, the STM tip would promote tunneling to/from 
the evanescent tail of dangling bonds. Therefore, cut¬ 
ting through the Si-Si bonds (at halfway between atomic 
planes) should provide a better description of the spatial 
dependence of the tunneling current in comparison, for 
instance, to cutting through atomic planes. As we shall 
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FIG. 2. Top: STM real space images of different donors 
showing the butterfly B (a) and caterpillar C (b) shapes. 
Bottom: images for the charge density from the wavefunc- 
tion of donors at (c) at z = 5.875 asi and (d) 
z = 5.625 asi. The color scale in hgures (c) and (d) ranges 
from 0 to 2 X 10“^ (ug.^). The arrows indicate the surface 
dimer direction. 

see, by comparison with experimental images in Sec. |ml 
this is precisely the case. 


III. REAL SPACE RESULTS 

Measurements were performed as described in Ref. [4] . 
Briefly, we directly measured the electronic ground 
state of shallow donors buried several nm underneath 
hydrogen-terminated ( 100 ) surfaces, by scanning tunnel¬ 
ing spectroscopy (STS) in ultra high vacuum (UHV). 
Samples with either P or As donors were prepared and 
measured. The As-doped samples contained donors lo¬ 
cated at random depths [4]. Samples with P donors were 
fabricated by submonolayer PH 3 dosing and well-defined 
encapsulation by epitaxial silicon, using a procedure sim¬ 
ilar to the one in Ref. m, but with an n-type substrate to 
promote elastic resonant transport HE]. Experiments 
on the donors were carried out with atomic resolution in 
real space, in the single-electron transport regime with 
both the tip and a heavily-doped region of the sample, 
at liquid Helium temperatures, acting as transport reser¬ 
voirs [ 12 ]. Both the As and P donors were determined 
to be in a lightly doped ^ 10 nm thick layer of silicon[4], 
with nominally one donor per 30 nm x 30 nm surface 
area. Analysis of dl/dU lineshapes, where U is the sam¬ 
ple bias and / is the tunneling current, revealed that the 
donors were typically coupled weakly to a buried reser¬ 
voir, relative to the thermal energy of excitations in the 


source and drain [TT] {ksT = 0.36 meV). 

To obtain high-resolution images of the neutral donor 
(D^) orbitals we adopted an unconventional scheme em¬ 
ploying an open-loop current measurement in the gap. 
The first scan maps the topography of the hydrogen- 
terminated (100) surface at U = —1.45 V, revealing 
dimerization along [110] directions. A second scan is 
made, following the topography of the first scan, but at 
a bias U such that only the lowest energy state is in 
the bias window {U 2 = —0.8 V). This scheme allows us to 
routinely image with very high resolution the donor or¬ 
bitals, whose surface electronic density can occupy more 
than 10 nm x 10 nm, with no cross-talk from the bulk 
states, or the two-electron (D“) localized state which is 
found at U ~ — 1.1 V sample bias. HE] 

The top panels of Fig. [^represent the two typical STM 
images for donors as obtained in our second scan mea¬ 
surements, with the dimer direction from the first scan 
indicated by an arrow. Both are twofold symmetric. Be¬ 
cause of their general shapes they are named here butter- 
fly (B) [Fig.^a)] and caterpillar (C) [Fig.j^b)]. The but¬ 
terfly is characterized by a nodal line in the low-frequency 
probability density that is always perpendicular to the 
dimer rows direction, while in contrast, the caterpillar is 
characterized by an antinodal line perpendicular to the 
dimer. 

We now contrast the experiments with results from the 
model described in the previous section. Figs. (a) to 
(d) show a sequence of inequivalent plane images calcu¬ 
lated from Eq. (j^ with 6 = 0, ie., planes with 

j = 0,..., 3. These correspond to each of the four in¬ 
equivalent atomic planes at d > 4, while panels (e) to (h) 
show the four inequivalent interstitial {5 = 1 / 8 ) cuts 
with j = 0,..., 3. Panel (a) is very similar to the well es¬ 
tablished behavior of the charge density cut for the ( 001 ) 
atomic plane containing the donor (see, for example. Fig. 
2(d) in Ref. 

Panels (a) and (c) in Fig. ^correspond to atomic planes 
belonging to the same sublattice as the donor. The cut 
images display fourfold symmetry (cross-like) around the 
impurity projected position. The atomic planes in the 
sublattice not containing the donor lead to lower sym¬ 
metry patterns (twofold, checkerboard-like). None of the 
atomic planes reproduce the STM images in the top pan¬ 
els of Fig. 

In contrast, calculated charge distributions at intersti¬ 
tial planes reveal very clear similarity to experimentally 
observed patterns (see Fig.|^. The experimentally mea¬ 
sured images in the upper panels compare very well with 
the theoretically generated images for interstitial planes 
and in the lower panels. These are represen¬ 
tative of B and C patterns, respectively. In general, the 
butterfly shape is found for interstitial cuts for j = 0 

or 3, z.e., just above or below an plane. The cater¬ 
pillar shape corresponds to images at for j = 1 or 2 . 
The comparison between experimental data and theoreti- 
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cal calculations is compelling, including the orientation of 
the charge distribution with respect to the surface dimer 
rows (determined by the direction of the dangling bonds, 
but not explicitly accounted for in theory). 

The calculated images corresponding to n = 5 and 
n = 4 /-planes (not shown) also manifest close similar¬ 
ities, but only on the overall shape (B or C). There are 
differences between the rapid oscillatory patterns, as ex¬ 
pected since the complete image is not periodic along 2 ;. 

We should stress that deviations from s-like charge 
densities have been observed for holes bound to Mn 
dopants in GaAs crystals in the form of a localized 
bowtie-shaped density. [14] This is a consequence of or¬ 
bital hybridization orbital hybridization in the KL impu¬ 
rity model. Such an effect is generic to the valence band; 
the signature of d-orbital hybridization can also be seen 
for boron acceptors on the [100] surface of silicon [T5]. We 
show here that the origin of the unusual low symmetry of 
the images in Ref. |4] is the true interference of valleys, 
not an edge effect. 


IV. UNDERSTANDING THE OBSERVED 

FEATURES: WAVEFUNCTION DISSECTION 

An intriguing result, related to the anisotropy of 
the donor bound state images, is the cyclic se¬ 
quence of butterfly (B) and caterpillar (C) shapes 
as the donor’s depth is varied [see Fig. i (e-h)]: 
...B(I)C(l)C(i)B(l)B(l),C(l)C(i)B(l)..., where S(i) 
means S-shape with axis aligned with [liO]. In words, 
B and C shapes alternate every two interstitial cuts and 
the symmetry axis changes from [110] to [110] for succes¬ 
sive images of the same shape. 

In order to investigate which element of the wavefunc- 
tion leads to this behavior, we compare the complete KL 
model with less accurate approximations. Fig. shows 
cuts across an interstitial B-type plane, /g^^, compar¬ 
ing the following cases: (a) The complete KL expression, 
Eqs. 0 and 0 , (b) isotropic effective mass, simulated 
by taking a = b = 0.815 nm in Eq. Q, (c) = 1, thus 

disregarding the periodic parts of the Bloch functions, or 
(d) eliminating the plane-wave parts in Eq. Q, as 

if k/^ = 0 for all /i. 

Anisotropic charge distributions are obtained in all 
cases except for the = 1 case, where a fourfold sym¬ 
metric image results. We conclude that the reduced sym¬ 
metry of interatomic planes on a diamond structure de¬ 
termines the observed fingerprint of anisotropy, since it 
is the periodic u that encodes these features. The same 
qualitative result is found in C-type images. In other 
words, the twofold symmetry of the images comes from 
the lattice. 

It is possible to understand qualitatively the STM fig¬ 
ures by simply analysing the geometric structure of bulk 
Si, see Eig. The B and C shapes result from the in¬ 
terference between the six valley states, imposed by the 
Ai ground state (see Eq. [^. Precursor shapes are clearly 



FIG. 3. Tomography of the KL wavefunction. Left: Cuts 
through atomic planes (a) (b) (c) ^ind (d) 

Right: Interstitial cuts (e) /q^\ (f) /j^^, (g) /j ^2 
(h) 43 / 4 . The color code represents |T( 2 : = zq)\^ in units of 

10~^(ag-^). The charge distributions calculated at the inter¬ 
stitial planes compare very well with the STM images, see 
Sec, mil for details. 


identified if we analyse an Ai-symmetric superposition 
v|/CBE(r) = T^e*'^--%(r). (4) 

This is illustrated in Fig. [^d), where we plot the elec¬ 
tronic density cuts |Tcbe(^, ^, ^/) P at interstitial planes 
for heights z/ = 6 + 1/8 and z/ = 6 + 1/2 + 1/8 

(^V2)- 

The 90-degrees rotation of the G and B symmetry axis 
in successive interstitial planes reflects alternation in the 
orientation of the zigzag pattern of the Si-Si bonds they 
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FIG. 4. Comparison between approximations for the charge 
distributions at plane z.e., at z = 5.125asi. (a) Full wave- 
function as described in Eq.Hwith anisotropic mass, (b) Full 
wavefunction with spherical (isotropic) mass, (c) Wavefunc- 
tion with trivial periodic part = 1. (d) Wavefunction for 
kfj^ = 0 for all fi. The twofold symmetry is maintained as far 
as the periodic part is kept. Equivalent results are ob¬ 
tained for the caterpillar shape. The color scale ranges from 
0 to 2 X 10“^ (ug.^). 


EIG. 5. Top: STM fourier space images of different donors 
showing the butterfly B (a) and caterpillar C (b) shapes. The 
four corners of the outermost dashed square boundary are 
reciprocal lattice vectors 27r(±l, ±l)/asi. Bottom: images for 
the theoretical effective mass wavefunction of donors located 
at z = 0 as seen from the interstitial planes (c) at z = 
5.125asi and (d) at = 8.375asi. 


cross [Fig. [^b)]. Interestingly, the orientation of the 
zigzag pattern of the Si-Si dangling bonds in the unrecon¬ 
structed surface plane also determines the orientation of 
the dimer rows in the reconstructed surface. Therefore, 
the orientation of surface dimerization and anisotropy 
of the donor STM image are completely correlated (al¬ 
though not related by causality - they share a common 
cause). In fact, from Figs.j^b) and (d) we note that the 
B or C symmetry axis orientation is always perpendicular 
to the zigzag pattern of the Si-Si bonds, which coincide 
with the dimer direction. 

As a final remark, we address the fast oscillatory fea¬ 
tures in the charge distribution images. Both the plane- 
wave and the periodic parts contribute to charge oscilla¬ 
tions as observed in Fig. The plane-wave components 
are the least intuitive, leading to oscillations that are in¬ 
commensurate with the lattice, as discussed in Ref. [4]. 
The envelopes simply confine the interference pattern 
around the donor position and do not contribute to the 
oscillations. 

Among the possible crystal orientations, [001] suffers 
from the least stringent oscillations of all directions, as 
pointed out in numerous works HU IT], resulting in the 
well determined cyclic sequence for (001) cuts in the 
donor charge distribution, following the lattice structure 
symmetry. This is related to the six phase factors 
for this particular orientation of the lattice. In a general 
direction, the periodicity is lost, see Appendix [A| 


V. UNDERSTANDING THE OBSERVED 
FEATURES IN THE FOURIER SPACE 


In this section, we use KL equations Q-Q to con¬ 
struct an analytic model describing quantum interference 
processes in the donor probability density |ff^(r)p. These 
interference processes ultimately originate from linear su¬ 
perpositions of lattice-incommensurate valley wavevec- 
tors and lattice commensurate reciprocal lattice vec¬ 
tors G. The model provides additional insight into the 
sequence B(1)C(1)C(I)B(1) in the interstitial planes 
for j = 0,1, 2, 3, and explains the origin of the sequence’s 
lattice periodicity. Moreover, the model forms the basis 
for the analysis in Section [VT| where we extract quantita¬ 
tive information about ko and h/a from the donor bound 
states. 

Writing u^{y) = exp(iG • r), we find that 

|ff^(r)p can be written as, 

|^(r)|2= ^ 


that is, slowly varying envelopes modulating oscillatory 
functions with spatial frequencies q = — k^+k^/ — G+G'. 
Here the summation is over —G and —k^ in ff^*(r) and 
G' and k' in ff^(r). Enumerating k^, k^/ and — G + G' 
in Table m we identify the associated two-dimensional 
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TABLE L Mapping of plane-wave components in T(r) and 
T*(r), to two-dimensional spatial frequencies ^in pd{x,y). 


class 

K 


-G + G' 
(27r/asi) 

q = iqx,qy) 

(27r/asi) 

1 


k^ 

(0,0) 

(0,0) 

1 



(0,0) 

(0,0) 

1 

±k^2 


(0,0) 

(0,0) 

2 

Tk/LiCC 

-\-^yz 

(0,0) 

qex± = (±0.85,0) 

2 


-\-^yz 

(0,0) 

qey± — (0, ±0.85) 

2 


-^yz 

(0,0) 

(0,0) 

3 

iky^x 

Tk/icc 

(±2,0) 

qpx± = (±0.3, 0) 

3 

^^yy 

F^yy 

(0,±2) 

<fpy± = (0,±0.3) 

3 

ikya,cc 

^^yy 

(±l.=Fl) 

qk± =±(0.15,-0.15) 

3 

iky^a, 

F^yy 

(±1,±1) 

qb± = ±(0.15,0.15) 


spatial frequencies q = (qx^Qy)- 

The main insight of Table |I] is to classify the oscil¬ 
latory frequencies in Equation (§, and to identify the 
inter-valley scattering processes producing special spa¬ 
tial frequencies q. The first class in Table [T| is the set 
of trivial non-interfering terms with q = 0 which oc¬ 
cur for + G' = + G. The second class is the 

set of interference terms whose oscillatory frequency q 
does not have lattice periodicity asi in the z direction, 
due to the presence of ±z valleys. Spatial frequencies 
qex± = ±27r(0.85,0)/asi and qey± = ±27r(0, 0.85)/asi be¬ 
long to this category, and originate from cross-terms of 
Tx and Ty valleys, respectively, with z and — z valleys. 
Since G = G', the terms listed in class 2 in Table |T] corre¬ 
spond to interference due to scattering by the valley-orbit 
potential of the donor. 

A third class in Table [I] is the set of interference terms 
involving the x and y valleys and G' 7^ G. These 
terms correspond to interference due to scattering be¬ 
tween valleys produced by the sharp central-cell poten¬ 
tial of the donor, and scattering by the lattice peri¬ 
odic potential. Since these terms do not involve the 
z valleys, their oscillatory frequencies q have the lat¬ 
tice period asi in the z direction. Spatial frequencies 
qpx± = 27r(±0.3,0)/asi and qpy± = 27r(0, ±0.3)/asi be¬ 
long to this category, which are cross-terms of ±x and 
Tx valleys, and cross-terms of ±y and =Fy valleys, re¬ 
spectively. Similarly, q]^± = ±27r(0.15, —0.15)/asi and 
qi)± = ±27r(0.15, 0.15)/asi belong to this category, and 
are cross terms of x and y valleys. 

We now turn to the real parts of the Fourier trans¬ 
forms of the data in Fig. [^a) and Fig. [^b), shown in 
Fig. [^a) and Fig. ib) respectively. As indicated in 
the colour scale, red-yellow (blue-cyan) denotes positive 
(negative) values. In both Fig. [^a) and Fig.j^b) we ob¬ 
serve ellipse-like features in the vicinity of qex± and qey± 
and more complex features in the vicinity of g = 0, as 
reported in Ref. [4]. The presentation of Fourier trans¬ 
forms herein differs compared to Ref. [4] only by showing 


the real parts, rather than the absolute value. 

Examining a large set of measured donors we find that 
the B and C symmetry charge densities differ in Fourier 
space along the direction gx (Fig- that we define as 
perpendicular to the dimer’s rows primary spatial fre¬ 
quency. For example, the B(l) pattern in Fig.a) is neg¬ 
ative at the g 5 + = 27r(0.15, 0.15)/asi side peak, positive 
at g = 0, and negative at the g^- = 27r(—0.15, — 0.15)/asi 
side peak. In contrast, the C(l) pattern [Fig. |^b)] 
is positive at the g/c+ = 27r(0.15, — 0.15)/asi side peak, 
positive at g = 0 , and again positive at the qk- = 
27r(—0.15, 0.15)/asi side peak. Features within gray 
dashed circles centred at 7r(±l, ±l)/asi are due to the 
2x1 reconstruction. [4] 

For comparison with theory we show the real part 
of the Fourier transforms for planes and in 
Fig. [^c) and Fig. [^d) respectively, which are B(I) and 
C(l) patterns. First we note that the features observed 
at qex± 9 .nd qey± in the measurements are reproduced, 
while reconstruction-related features are absent, because 
of assumption (in) in Section]^ Focusing on the region 
l^ccU^I < 7r/asi, the g^^ side peaks are negative for the 
B(I) pattern, and the qk± side peaks are positive for the 
C(l) pattern. Consequently, the interstitial plane repro¬ 
duces the main Fourier-space features classifying the B 
and C patterns in the measurements. 

In rows 1 and 2 of Fig. [^a) we present in order the se¬ 
quence of Fourier transforms alternating between atomic 
and interstitial for n = 5. First, we note that the in¬ 
equivalence between [ 110 ] and [ 110 ] directions in the in¬ 
terstitial planes is reflected by alternation of the sign and 
amplitude of peaks at g^^ and qk±. As expected from q 
in Table [Tj the peaks centred at g^± and qk± are periodic 
with the lattice constant asi in the z direction. Peaks 
at qpx± and qpy± are absent in the interstitial plane, but 
represented in the atomic plane, and as expected from q 
in Table [I| they vary twice as fast with 2 ; compared with 
qb± and qk±. 

VI. CONDUCTION BAND MINIMA AND 
ENVELOPE ANISOTROPY 

The anisotropy 6 /a of the envelope functions and the 
wavevector |k^| = ^0 of the conduction band minima in 
KL equations Q and <§ can be determined from the 
spatial frequency and anisotropy of the ellipse-shaped 
features found in the Fourier space data at qex± = ±x/co 
and qey± = ±yA^o- Our analysis relies on the isolation of 
these features in Fourier space, which are highlighted for 
a P donor in Figure [^a), and the specific origin of this 
interference term as identified in Table [H 

Owing to localization of the donor bound states, spa¬ 
tial oscillations at wavevectors q = —k^ + k^/ — G + G' 
in Table [I] and Equation modulate slowly varying en¬ 
velopes r)Fp{r). Focusing on the rele¬ 

vant k^, G, and G' in Table IS we determine from 
Equation ^ that contributions to pT(r)p at frequencies 
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FIG. 6. Interference processes for the atomic and interstitial planes in the range \qx\ ^ 27r/asi and \qy\ < 27r/asi. 



FIG. 7. Top: STM Fourier space image of a caterpillar (G) 
shape (a). The four corners of the outermost dashed bound¬ 
ary are reciprocal lattice vectors 27r(±l, ±l)/asi- Bottom: 
Fourier filtered data isolating the x (b) and y (c) ellipse fea¬ 
tures (top-left inset), real space data shifted back to the origin 
in Fourier space, and best ht to bulk-truncation model from 
the main text. 


Qex± = ix/co and qey± = iy^o are of the form 

nex{r) = Cx cos{koz) cos{kox)Fx{T)Fz{r), and (6) 
ney{r) = Cy cos{koz) cos{koy)Fy{T)Fz{r). (7) 

For fixed z = neyir) is a product of envelope functions 
modulated by oscillations cos(k^-r), with a constant pre¬ 
factor Cycos{zod), where depends on the Ak,G- 
Our model for the conduction band minima, described 
in Appendix!^ relies only on the definite parity of the 
donor envelope functions, and allows for the extraction 


of ko separately for the x and y direction. Fitting the 
data in Figure [^a) around q = x/cq and q = y/co sep¬ 
arately, we obtain essentially identical extrema k^ = 
(0.83 ± 0.02)(27r/asi) for the x and y conduction band 
valleys, respectively. The errors are dominated by inac¬ 
curacy in the determination of q relative to lattice fre¬ 
quencies 27r/asi, since the latter are blurred in topogra¬ 
phy due to finite sampling and instrumental errors, as 
discussed in Appendix [B| 

The observed anisotropy of the Fourier representa¬ 
tion of necc(r) and ney(r) in Fig. [^a) results from the 
anisotropy h/a of Fx{x^y,d) and Fy{x,y,d) respectively, 
since Fz{x^y^d) is isotropic in the x — y plane. To de¬ 
termine 6/a, we first isolated the ellipse-like features in 
Figure ^a) corresponding to nex{^) and ney(r) in our 
model by applying a filter in Fourier space passing the 
relevant data inside the green boundaries. Results of the 
filter for the x and y valleys are shown in the upper-left 
insets of Figure [^b) and Figure [^c), respectively. Now 
isolated, these features can be directly fit to model ex¬ 
pressions © and 0- 

Independent of the details of the Fourier filter em¬ 
ployed to isolate the features, we obtain envelope 
anisotropies h/a = 0.54 ± 0.02 and h/a = 0.53 ± 0.02 for 
the X and y valleys, respectively. The real space data are 
shown in Figure [^b) and Figure [^c) after shifting back 
to the origin in Fourier space using ko from above. The 
anisotropy of this data reflects the anisotropy of the en¬ 
velopes Fx{v) and Fyir). The lower right insets of Figure 
l^a) and CTb) respectively shows the model calculation 
for the least-squares parameters, demonstrating excellent 
agreement with the data. 

The results obtained for the donor in Figure [^a) are 
representative of other donors. Across a total of six 
donors (four As donors and two P donors), the val¬ 
ues obtained for ko are very reproducible; best fit val¬ 
ues for ko in the range [0.81, 0.84] (27r/asi) are obtained, 
with similar confidence intervals between 0.02(27r/asi) 
and 0.03(27r/asi). The mean value of all measurements 
is ko = (0.82 ± 0.03)(27r/asi). Across the same set of six 
donors, we obtain h/a = 0.49±0.03 as the mean value for 
both the X and y valleys. Nine of the twelve best fit values 
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for hja fall in the range [0.47, 0.58], though we also found 
three anomalously smaller values of ^ 0.40 ± 0.03. The 
larger spread of 5/a values could be caused by disorder 
induced by the presence of randomly located dangling 
bonds on the surface. We observe that these dangling 
bonds, whose position can be identified with atomic res¬ 
olution, spatially modulate the density of states of the 
conduction and valence bands on length scales similar to 
the donor’s effective Bohr radius. 


VII. DISCUSSIONS AND CONCLUSION 

The notable resemblance between the conduction 
splotches and a cut of the KL wavefunction charge den¬ 
sity entails two results. Firstly, that the KL theory is not 
limited to the estimation of the energetics of an electron 
bound to a dopant. It may be used as a good model of the 
electron wavefunction as well. Secondly, that the com¬ 
plex profile of the STM image is a surprisingly accurate 
depiction of the bulk wavefunction, and the surface does 
not disturb the electronic wavefunction beyond recogni¬ 
tion. 

From the first result, we may extrapolate important 
implications to the design of quantum devices. Theo¬ 
retical proposals based on the KL model of the donor 
wavefunction should be reasonably accurate, as long as 
the fast oscillations of the wavefunction are consistently 
taken into account liiniiiH]. Moreover, the resilience of 
the valley coherent interference and overall shape of the 
wavefunction reiterates that the donor-bound electron is 
robustly shaped by the donor potential and valley-orbit 
coupling, and disturbances like the passivated interface 
and the STM tip do not alter significantly the dopant 
wavefunction. 

From the second result we conclude that the reduced 
symmetry of the STM images is not an artifact of the 
experiment, but reflects the true nature of the bulk donor 
wavefunction. Measurements made at the surface are 
good estimates of the bulk behavior. We may use that 
to appoint experimental values to the wavevector at the 
conduction band minima ko = \k^\ and the anisotropy 
hja. 

While early reported measurements for ko vary from 
0.77 ^'Kjasi [E] to 0.85 27r/a5i[20], we obtain ko = 
(0.82 ± 0.03)(27r/asi), which is independent of the model 
adopted for the envelope function of the dopant ground 
state. 

Adopting a KL envelope, we obtain a ratio b/a = 
0.49 ± 0.03 that fits the experimental data. Most fre¬ 
quently, we measure ratios near 0.52, in excellent agree¬ 
ment with the theoretical values obtained within a cen¬ 
tral cell corrected model (0.53 for both As and P) and is 
slightly lower than the ratio obtained by KL without cen¬ 
tral cell (0.58 using modern values of the effective mass 
and dielectric constant). It is expected that b/a changes 
as the wavefunction is probed at diferent distances from 
the impurity center. Near the center, the Schrodinger 


equation is dominated by the nearly spherically symmet¬ 
ric potential energy, pushing the ratio b/a towards 1. Far 
from the center, however, the kinetic energy dominates 
at large, and one expects b/a ^ ^/Ta±Jrr^ = 0.46. This 
could explain the fits from measurements being system¬ 
atically lower than the theoretical estimates. 

Ideally, the anisotropy would be estimated both the¬ 
oretically and experimentally as a function of the dis¬ 
tance from the nucleus. Unfortunately, a trial variational 
wavefunction with too many fitting parameters leads to 
unreliable results, and the lack of a definite estimate of 
donor depth impairs this connection from the experimen¬ 
tal point of view. A rigorous description should allow for 
the anisotropy to be a function of the distance from the 
nucleus. Unfortunately, within the variational scheme, 
this involves increasing the number of parameters in the 
trial envelope functions, with questionable results. At 
the current stage of Si quantum technologies [3], most 
applications do not urgently require such refinement. 

The KL theory is shown here to be accurate and appli¬ 
cable to new quantum technologies - a remarkable feat 
for a 60-year-old model. 
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Appendix A: Perspectives for sub-(lll) surface 
donors 

Experimental data presented herein focuses on the 
H/Si(100)2xl surface, prepared by UHV annealing to 
produce large defect free terraces whose dangling bonds 
are saturated by UHV exposure to atomic hydrogen. [21] 
While preparation of hydrogen terminated Si(llO) sur¬ 
faces is more difficult, [22| the Si(lll) surface presents 
another immediate possibility since large defect-free ter¬ 
races of the well-known H/Si(111)7x7 surface can be ob¬ 
tained by aqueous treatment in hydrofluoric acid. [23| 
The 1 X 1 reconstruction of H/Si(lll) can also be ob¬ 
tained by wet chemical treatment, [24| or in UHV by 
thermal treatment of the H/Si(111)7x7 surface under 
flux of atomic hydrogen. [25| 

We present in this appendix images to be expected for 
STM of donors under surfaces with (111) orientation. In 
analogy with the results for (001) surfaces, one might 
expect a cycle of 3 inequivalent charge distributions at 
successive (111) interstitial cuts, which is not obtained 
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FIG. Al. Charge densities calculated at two structurally 
equivalent interstitial [111] plane at a distance 4.125 x 
from the donor. The color scale ranges from 0 to 5 x 10 
on both images. Note significant differences at the center and 
close similarity toward the edges due to the severe valley in- 
tereference along the (111) direction. 


FIG. A2. Butterfly images of the charge distribution from 
donors at two equivalent (001) cuts but at different depths 
(a) and (b) The more distant donor (b) spreads 

over a larger area of the surface. The color scale ranges from 
0 to (a) 4 X 10“^(ag.^) and (b) 10“® x 


here. We observe all cuts with an overall triangular shape 
at large scales and a central part with no periodicity or 
similarity among different cuts for all donor depths we 
have examined. 

We attribute the absence of periodicity in the central 
region of the figures to the interference effects among 
the incommensurate plane waves at and the lattice- 
commensurate plane waves expansion of the periodic 
functions u^. Interference effects are known to be 
stronger along (111) than (001) directions [16]. 

Appendix B: Uncertainty in ko 

Real space confinement (5r ^ 3 nm of surface charge 
density of the bound state blurs the Fourier space 
peaks by an amount Sq ^ 27rlSr ^ 0.2(27r/asi). 
Nevertheless, ko can be determined in Fourier 
space with accuracy better than 6q, as follows, 
because of the inversion symmetry of the donor en¬ 
velope functions. Expanding Equation around 


q = x/co in Eourier space, we obtain nex±{(lx^(ly^ z) = 
Y.n,m^n,m{qx ± koYq^/n\m\, where Mn,m = 
\[d'^+^Fx:,{q^,qy,z)/dq^dq^]\ij=o, and Fxz{qx,qy,z) = 
f f z)Fz(x,y, z). Values for 

^ f f dxdy(-ix)^(-iy)"^Fx(x,y,z)Fz(x,y,z) 
are non-zero when n and m are both even, when the 
integrand has even parity in x and y coordinates. To 
lowest non-zero order in q, we have riexiqx^Qy^ = 
Mofi - \M 2 fi\{qx ± koY/2 - |Mo, 2 |g ^/2 in the vicinity 
of g = x/cq. Eollowing the same approach, we obtain 

ney{qx,qv,z) = Mo,o “ 1-^2,o|(% ± fco)^/2 “ |Mo,2|g^/2 

in the vicinity of g = yko. In other words, define the 
extrema in the Eourier space distributions of Equations 

and [^. 


Appendix C: Estimating the donor position 

We investigate to what extent the image variety pre¬ 
sented here may provide information on the (x, z) po- 
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sition of isolated substitutional donors in the Si lattice. 
If feasible, this would guide sorting favourable samples to 
perform specific spin qubit operations. For example, in 
a Si-donor-based quantum computer inspired by Kane’s 
original idea, [26] spin qubit operations are known to be 
highly sensitive to the donor position m- In particu¬ 
lar, two-qubit exchange gates are anticipated to require 
precise positioning of the interacting donor pairs. 

All theoretical images have a center of inversion sym¬ 
metry, which readily identifies the [100] and [010] (or x 
and y) coordinates of the donor. This is not the case 
for the STM images, where patterns may or may not 
have a center of inversion symmetry (see Fig. |^a) and 
(b) respectively). This aspect is related to the relative 
(x, y) positioning of the donor with respect to the surface 
dimer rows, which may break one of the mirror symme¬ 
tries (tranversal to the caterpillar figure) of the STM im¬ 
ages. Properly taking into account this effect could lead 
to the donor (x, y) position. 

Laboratory fabrication of buried dopants in Si with 
atomically precise z-locations relative to the surface is a 
complicated issue at this early stage. Uncertainties may 
be introduced by dopant diffusion during annealing. A 
posteriori analysis of the shift of the conduction band 
edge can give a clue on the dopant depth with an error 
of the order of the lattice parameter. [4] 

If the donor depth is determined by least-squares fit¬ 
ting of the band edge potential to an interval [d — asi, d + 
asi] containing 8 planes, with 95% confidence, then the 
categorization into a B or C-type image reduces the num¬ 
ber of possible planes to 4. 

Note that all features highlighted here are specific to 
(001) surfaces. In principle, the (111) surfaces discussed 
in Appendix form kaleidoscope-like figures that could 
be used to distinguish the donor depth uniquely. 

Although all the results presented here correspond to 
a single donor, we may try to extract information about 
the relative position of two donors. A simple check for the 
same class of images is to compare the range of the charge 
distribution probed by STM. More distant donors corre¬ 
spond to more spread images, as illustrated in Fig. |A2| 


Identical images correspond theoretically to donors ly¬ 
ing on the same xy plane. However this is not a strict 
experimental requirement - two donors with the same 
^ coordinates may appear differently at the surface due 
to the Si dimer rows relative position with respect to 
each of them. Therefore, identification of donor pairs 
at the same depth below the surface requires careful con¬ 
sideration of the structural peculiarities of Si consistently 
combined with the reconstructed surface - each donor’s 
position relative to the surface dimers must be included 
in the analysis of the STM images. This feature may 
be included in our model, but the number of possible 
combinations do not fit into simple rules and should be 
analysed systematically in each case. 

Appendix D: Oscillatory behavior and donor-based 
qubits 

It is possible to harness spins for quantum computa¬ 
tion if we are able to control the spin-spin exchange cou¬ 
pling on demand. This task is difficult due to the sensi¬ 
tivity of the exchange coupling J to the donors relative 
positioning miiniiTH]. Exchange coupling oscillations 
arise from the interference between the plane wave parts 
of two donors wavefunctions at a relative position R. 
Each donor establishes a pinning point for the 6 plane 
waves, leading to factors of the form exp[i(k^ — * R] 

in the exchange coupling. As a result rapidly oscillating 
coupling J(R) is obtained along general directions R. 
There is however one favourable situation: donors kept 
at R strictly along one of the (100) directions lead to a 
smoother behavior for J(R). 

The fast oscillatory patterns of charge in our theoret¬ 
ical or experimental images are related not only to the 
plane wave part of the Bloch functions, but also show 
interferences coming from the periodic part of the Bloch 
functions. It is clear that the oscillatory behavior of the 
exchange coupling J is not the same as the oscillatory 
pattern of the electronic wavefunction. The exchange os¬ 
cillations are only due to the valley interference, since the 
periodic u^(v) always interfere constructively. 
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